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Abstract 

This paper is concerned with the numerical solution of a class of variational inequalities of 
the second kind, involving the p-Laplacian operator. This kind of problems arise, for instance, 
in the mathematical modelling of non-Newtonian fluids. We study these problems by using a 
regularization approach, based on a Huber smoothing process. Well posedness of the regularized 
problems is proved, and convergence of the regularized solutions to the solution of the original 
problem is verified. We propose a preconditioned descent method for the numerical solution of 
these problems and analyze the convergence of this method in function spaces. The existence 
of admissible descent directions is established by variational methods and admissible steps are 
obtained by a backtracking algorithm which approximates the objective functional by polynomial 
models. Finally, several numerical experiments are carried out to show the efficiency of the 
methodology here introduced. 
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1 Introduction 

Variational inequalities (Vis) provide a versatile background for the analysis and modelling of phys¬ 
ical phenomena which involve free boundary problems. This kind of problems include, for instance, 
contact of rigid bodies, flow of electro- and magneto-rheological fluids and flow of viscoplastic ma¬ 
terials, among others (see H1H2II31IEU25I). The wide range of applications of this kind of free 
boundary problems make the analysis and the numerical simulation of their associated Vis a quite 
interesting and challenging held of research. 

On the other hand, the p-Laplacian operator has been widely analysed as a model case for 
quasilinear and degenerate elliptic equations (see PES]). Regarding the p-Laplacian problem, 
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several analytical results, concerning existence and multiplicity of solutions, have been obtained in, 
e.g., [381139| . Further, the numerical analysis of the p-Laplacian problem has been a productive field 
of research. Mainly, the finite element approximation has been broadly studied in the literature in, 
for example, mmm and the references therein. The numerical realisation of the p-Laplacian has 
been carried out by the Augmented Lagrangian method m and, recently, by using optimization 
and variational techniques |25j and multigrid algorithms |3j. 

In spite of the fact that the p-Laplacian operator is a extensively studied field, scarce work can 
be found in the numerical analysis of variational inequalities involving this differential operator. 
The obstacle problem with associated operator of the p-Laplacian type has been analyzed in [ 32] . 
There, the authors analize a finite element approximation of the problem and provide error estimates 
for such approximation. In [28], the obstacle problem in the context of a glaceology application 
has been studied. The authors consider a variational inequality of the first kind, involving the 
p-Laplacian operator. They analyze the existence, uniqueness and regularity of solutions, and 
propose a finite element approximation of the problem. Elliptic and parabolic quasi-variational 
inequalities involving quasilinear operators are considered in [22, '23]. In these papers, the authors 
propose and study a semismooth Newton approach for the numerical solution of these problems, 
and provide several theoretical results regarding existence and regularity of solutions. Finally, 
in [25], the authors propose a finite dimensional descent algorithm for the numerical solution of 
several differentiable problems, including the classical Dirichlet p-Laplacian problem and a class of 
variational inequalities wich combines the Laplacian and the p-Laplacian operators, for 1 < p < oo. 

In contrast with the previous contributions, in this paper we are concerned with variational 
inequalities of the second kind. The importance of this class of variational inequalities lies on the 
fact that they can be used to model the flow of a particular class of viscoplastic materials: the 
Herschel-Bulkley fluids. 

Herschel-Bulkley is a power-law model with plasticity. This model is used to simulate some 
materials whose behaviour depends on the flow index p. This constant measures the degree to 
which the fluid is shear-thinning (1 < p < 2) or shear-thickening (p > 2). The Herschel-Bulkley 
model can be seen as a generalization of the classical Bingham model, which is retrieved from the 
first one by taking p = 2 ([32, EH)- Depending on the value of the power index, this model can 
be used to simulate a wide range of materials from nail polish or whipped cream (shear-thinning 
fluids) to quicksand or silly putty (shear-thickening fluids) (see 0). Furthermore, the Herschel- 
Bulkley model has proved to be accurate in the modelling of blood, a known shear-thinning fluid 

m S3 si]. 

In consequence, the numerical resolution of Vis involving a p-Laplacian operator is an important 
research field. A classical approach to these problems is the Augmented Lagrangian method (see 
[26]), while, from our point of view, the application of optimization and variational techniques has 
not been explored enough in this context. Several optimization problems involving non differentiable 
functionals have been successfully analyzed by using this approach (see m- Further, the analysis 
of this kind of methods in function spaces is a challenging but promising research field. 

As stated before, in this paper we are concerned with a class of variational inequalities of 
the second kind involving the p-Laplacian operator and the L 1 -norm of the gradient. Our main 
intention is to develop an efficient algorithm for the numerical solution of these problems. The 
main challenge in this aim consists in designing a numerical strategy, which allows to obtain an 
accurate solution with a fast convergent method. In the context of numerical solution of Vis of the 
second kind, the local smoothing techniques, such as Huber regularization, have proved to be an 
effective way to achieve such a goal (see dans]). Therefore, we study the variational inequality 
as an equivalent minimization problem of a non differentiable functional, and we regularize the 
functional by a Huber procedure. Further, the convergence of the regularized solutions to the 
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original one is established. 

For the numerical solution of the Vis under study, we propose a preconditioned descent algo¬ 
rithm in function spaces. Several issues arise in this approach. Mainly, we need to discuss the 
existence of admissible search directions and admissible step sizes. Admissibility of step sizes de¬ 
pends on the line search strategy. Here, we propose a backtracking algorithm which approximates 
the objective functional by polynomial models. In this way, the algorithm provides admissible step 
sizes with low computational effort. 

The existence of admissible search directions is analyzed considering the two cases 1 < p < 2 
and p > 2, separately. In the case 1 < p < 2, we first discuss the properties of a suitable Hilbert 
space in which we will propose and analyse the algorithm. Next, we define the preconditioner and 
prove the existence of admissible search directions. This is achieved by discussing the conditions 
for the Zoutendijk condition to hold ([33, JOj). Finally, we state and prove a global convergence 
result for this algorithm. The case p > 2 poses analytical issues which prevent us from studying 
the algorithms in function spaces. In fact, it is not possible to prove existence of admissible search 
directions in the same function space in which the elliptic preconditioner is defined. We discuss in 
detail these issues and propose an alternative algorithm in a finite element space. Next, we state 
and prove a global convergence result for this algorithm as well. 

Though the descent algorithms are usually slow, in our case the design of suitable precondi¬ 
tioners and the use of an innovative line search algorithm help us to obtain a robust algorithm 
which only needs the solution of one linear system per iteration. Further, since the algorithms are 
proposed, at least in the 1 < p < 2 case, in function spaces, they are expected to exhibit mesh 
independence. 

The paper is organized as follows. In Section [2] we introduce and analyze the variational inequal¬ 
ity and its associated optimization problem. Next, by using Fenchels duality theory, a necessary 
condition is derived and, since the original problem is ill-posed, a family of regularized optimiza¬ 
tion problems is introduced and the convergence of the regularized solutions to the original one 
is proved. In Section [3] the numerical approach to these problems is studied. We propose pre¬ 
conditioned descent algorithms for the regularized optimization problems, considering separately 
the two cases 1 < p < 2 and p > 2. Particularly, we prove a global convergence result for all 
the algorithms constructed. Section [4] is devoted to the numerical experience. First we discuss 
the main issue regarding the implementation of our algorithms. Mainly, the discretization issues 
and the implementation of the line-search methods. Next, several numerical experiments, which 
illustrate the main features of the proposed approach, are carried out. Finally, in Section [5j we 
outline conclusions on this work and discuss some challenging issues that can be analyzed in future 
contributions. 

2 Problem Statement and Regularization 

Let us start this section by introducing some important notation. The scalar product in and the 
Euclidean norm are denoted by (•, •) and | • |, respectively. The duality pairing between a Banach 
space V and its dual V* is represented by (■,-)v*,v> and || • ||y stands for the norm of V. Given 
1 < p < oo, the conjugate exponent is denoted by p'. Further, the duality pairing between L p and 
L p spaces is denoted by (•, -) p \ p (see [5J Th. 4.11]). We use the classical notation for the Sobolev 
space Wq’ p (Q) and the notation VF _1,p, (H) for its dual space. Finally, we use the following bold 
notation L P (H) := L P (Q) X L P (H) for 1 < p < oo. 

This work is concerned with the numerical solution of the following class of variational inqualities 
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of the second kind: find u E VF 0 '’ p (0) such that 

f |Vu| p ~ 2 (Vu, Vw) dx + g f \Vv\dx — g f \Vu\dx> [ f[v — u)dx , Vv G W 0 1,!! (fi), 

Jsi Jtt Jo, 

where 1 < p < oo, g > 0 and / E L p ' (LI). 

It is well known that this variational inequality represents a necessary optimality condition for 
the following optimization problem of a non-smooth functional 


min J(u) := - f \\7u\ p dx + g f |Vu| dx — f fudx. 

P dn Jn Jn 


( 2 . 1 ) 


Therefore, we will focus on the numerical solution of (2.1), by using optimization and variational 
techniques. 

Theorem 2.1. Let 1 < p < oo. Then, problem (2.1) has a unique solution u E 1Tq’ p (11). 

Proof. Note that functional </(•) can be rewritten as 

J(u) = -\\u\\ p 1:P +g [ | Via | dx — [fudx. 

P ^0 Jn Jn 

Therefore, it is clear that </(•) is a continuous and strictly convex functional, which satisfies that 

lim J(u) = Too. 


This fact yields (see, for instance, m Ch. 2] and [3T( Ch. 1]) the existence of a unique solution 
u E TITq’ P (S 7) for the problem (2.1). 


□ 


□ 


2.1 A Multiplier Characterization 


In this section, we use the Fenchel’s duality theory to characterize the solution of problem (2.1) 
with a vectorial function, which acts as a multiplier. The aim of such a procedure is to obtain an 
optimality system which will be used to characterize the solutions of (2.1). 

For the sake of readability of the paper, let us briefly describe the main ideas in Fenchel’s theory. 
Let V and W be two Banach spaces with dual spaces V* and W* , respectively. Let A E C(V, W) be 
given and let A* E C(W *, V*) be its conjugate functional. Further, let T : V —>■ M and Q : W —> M 
be two given functionals. We are concerned with minimization problems in which the objective 
functional J : V — > M can be decomposed as 

J{u) := T(u) T (7(A u). 

In such a case, the problems we are interested in are given by 


inf{.F(u)T£(Au)}. 

u&V 


Next, it is known that the associated dual problem of (2.2) is given by (see PH pp. 60-61]) 

sup {—-F*(— A*q) — C?*(q)}. 
q ew* 


( 2 . 2 ) 


(2.3) 
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Here, F* : V* —> M and G* : W* —> M denote the convex conjugate functionals of F and Q, 
respectively, i.e,, 


.F*(-A*q) = sup{(-A*q, v) v *y - F(v)} and C/*(q) 
vev 


sup {(q, r) w * >w - C/(r)} . 
r£W* 


Now, let us suppose that the primal problem (2.2) has a unique solution u £ V and that both F 
and Q are convex and continuous. Then, [El Th. p. 59] and [13 Rem. 4.2, p. 60] imply that no 
duality gap occurs, i.e., 


+ <3(Ku)} = sup {—J"*(—A*q) - <T(q)}, 

U&V q £W* 


and, moreover, that the dual problem has at least one solution q £ W*. 

Finally, Fenchel’s duality theory allows us to characterize both the primal and dual solutions. 
Indeed, d3 p. 61] implies that u and q satisfy the following system of equations 


— A*q £ dF{u) 
q £ dQ(Vu), 


(2.4) 

(2.5) 


where dF(u) and dQ(Vu ) stand for the subdifferential of F at u and the subdifferential of dQ at 
Vrt, respectively. 

Let us turn our attention to problem (2.1). First, we define V := H / 0 1 ’ P (H), V* := W~ l ' p '{£!), 
W = L P (H), and we identify the dual space of L P (H) with L p, (f2) (see (4| Th. 4.11]). 

Next, we introduce the functionals T : H / 0 1,P (H) -> I as F(u) := | f n |Vtt| p dx + f n fudx and 
Q : L p (12) —> M as ^(q) := g |q| dx. It can be easily verified that these two functionals are convex, 
continuous and proper. We also introduce the linear operator A : Wq’ p —> L p (fl) by A u := Vu. 
Clearly, A £ C{Wq , L P (H)). Thanks to these definitions, it is clear that problem (2.1) satisfy all 
the requirements of Fenchel’s duality theory. Therefore, there exists at least one solution for the 
dual problem. Moreover, the solutions of primal and dual problems u £ Wq ,p (Q) and q £ L p, (H), 
respectively, satisfy the system (2.4)-(2.5). 

First, we study (2.4). In this case, since T is Gateaux differentiable, the subdifferential of T 

al F' (see 

W’- 1 -p,w 0 1, 


reduces to the Gateaux differential F' (see H3 Prop. 5.3, p. 23]). Therefore, ( |2.4[ ) implies that 
(—A*q, v) w -i, P W } ,p = (F'{u) , v) w _ 1 , P W i,p, \/v £ W 0 1,p (O), which is equivalent to 


-(q,Vu) p> = [ \Vu\ p ~ 2 (Vu,Vv)dx- [ fvdx, \/v £ W 0 1,P (H). 
isi Jvi 

Now, thanks to the Riesz’s representation theorem in L p spaces (see [SJ Th. 4.11]), there exists a 
unique w £ L p (Q) such that (q, \/v) p ^ p = f n (w,Vv) dx, which yields that 

— f (w ,Vv)dx= I |Vu| p_2 (V : u, Vv) dx — f fv dx = 0, Vu £ H / 0 1,P (H). 

Jq Jfi Jq 

Next, we analyze (2.5). In this case, since Q is not differentiable, ( |2.5[ ) implies that 

G(Vu) - Q( r) > (q, X7u-r) p f )P , Vr £ L P (H). 

By following similar argumentation as in m Sec. 2.1], we conclude that the last expression implies 
that 

(q, r) p / iP <g |r| dx,Vr £ L P (H) and g / |Vw| dx = (q, r) p / iP . 

Jtt J f2 
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( 2 . 1 ) 


Finally, thanks to m Lem. 2.1] and [121 pp. 85], we obtain the following optimality system for 

f |Vu| p_ 2 (Vu, X7v) dx + f (w, Vu)dx — f fv dx = 0, Vv £ Wq’ p (D), (2.6a) 

|w(x)| < g, a.e. in Q, ( 2 . 6 b) 


\7u(x) = 0 or 
VF(x) / 0 and w(x) = 

Definition 2.2. The active and inactive sets of the problem are defined by 

A := {x £ 12 : Vu(x) / 0} and Z := {x £ fl : Vu(x) = 0}, 

respectively. 

2.2 A Huber Regularization Procedure 


( 2 . 6 c) 


The non-differentiability of problem (2.1) can provoke instabilities in several numerical schemes, 


such as a primal-dual algorithm (see na). This issue can be appreciated in the fact that system 


(2.6) does not have a unique solution. Further, this lack of regularity prevents us from developing 


an algorithm based on optimization techniques, as proposed. A classical approach to this kind of 
problems is regularization. However, the question about what kind of regularization procedure is 
the most suitable is a hot topic (see mmm)- 

In this work, we propose a local regularization of Huber type. The big advantage of using such 
a procedure is that Huber regularization only changes locally the structure of the functional in 


(2.1), preserving most of the qualitative properties of functional J. 


Let us start by introducing, for 7 > 0, the function by 


by 


i/jJz) : = 


g\z |-fy if 7^1 >9 


if 7 |z| < g. 


(2.7) 


Note that Vb corresponds to a local regularization of the Euclidean norm. In Figure [l] it is possible 
to appreciate the effect of this regularization in dimension one. 


Next, by using the function we propose the following regularized version of problem (2.1) 

min J 7 (it) := - j \Vu\ p dx+ j ip^(Vu)dx— j fudx. (2-8) 

«eVK 0 1 ’ p (O) V Jn in Jn 

Theorem 2.3. Let 1 < p < 00 and 7 > 0. Then, problem (2.8) has a unique solution £ ^^’^(H). 


Proof. First, let us state that function is a convex function m • Therefore, the functional J 7 (u) 
has the same qualitative properties of functional J(u). Consequently, the result follows in the same 
way as in Theorem 2.1 □ □ 


We again propose the use of Fenchel’s duality theory to generate an optimality system for 
Actually, in this case we only need to replace the functional Q by the functional t/ 7 : L p (f2) — > M 
given by ^(p) = Jq Vb(p) dx, which is convex and continuous. Therefore, we can use the Fenchel’s 
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Figure 1: Huber regularization in dimension 1. 


theory to state that the dual problem has at least one solution q 7 G L p, (H), and, moreover, that 
Uy and q 7 satisfy the system (2.4)-(2.5). 

Since the functional T has not changed, in this case (2.4) reads as follows 


— (q 7 , Vr)p- p = / |Vu 7 | p- 2 (V'u 7 , Vv) dx- fvdx, Vv G Wq P (Q). 

J J Q 


Further, thanks to the Riesz’s representation theorem in L p spaces (see [H Th. 4.11]), there exists 
a unique w 7 G L p (H) such that (q 7 , Vv) p ' tP = f n ( w 7 , Vv) dx. This fact yields that 


— / (w y,Vv)dx= / \Vuy\ p 2 {Vuy,Vv)dx — / fvdx = 0, Vu G Wq’ p (Q). 

J J J 

On the other hand, the functional Gy is Gateaux differentiable. Therefore, in this case equation 

(q 7 > r )p',p (&y(Vrt 7 ) , r) p / )P , 


(2.5) is given by 


which is equivalent to 


(q 7 > 


n f n 


7 g - dx , Vr G I/(H). 


' y max(£f,7|Vrr 7 |) 

Finally, since w 7 is the unique Riesz representative of q 7 , we have that 

[ (w 7 , r) dx = [ 75 - (V?i 7 ,r) ^ Vr G L P (H). 

Jo Jn max(5,7|Vu 7 |) 


Summarizing, we have the following regularized optimality system for (2.8). 


|Vu 7 | p 2 (Vuy, Vv) dx + / (w y,Vv)dx— / fvdx = 0, Vu G VFq’ p (H). 

«/ Q J Q, 


w 7 (x) = 57 


VUy(x) 


max( 5 , 7 |Vfx 7 (x)|) 


a.e. in H and 7 > 0 . 


(2.9) 


(2.10a) 

(2.10b) 
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Definition 2.4. The regularized active and inactive sets are given by 


A~ := {x £ Q : 7 |Vu 7 (.t)| > g} and := {x £ Ll : 7 |Vu 7 (x)| < g}, 


respectively. 


Lemma 2.5. Let 1 < p < oo and 7 > 0. Then, the sequence of optima of (2.8) is bounded in 

Wo’W 

Proof. Let us start by noticing that 

- [ |V« 7 | p Xr — j fu 7 dx < J 7 (u) < J 7 (0) = 0. 

P Jn Jn 

Next, Holder and Poincare inequalities imply the existence of a positive constant C , which only 
depends on Ll and p, such that 

\\\ u iW v w i,v < c\\f\\ LP '\\Vuj LP . 


p 


Since p > 1, the last expression directly implies the result. □ □ 

Theorem 2.6. Let 1 < p < 00 . Then, the sequence {u 7 } C W^ p {fL) converges strongly in Wq P (D) 
to the solution u of problem ( 2 . 1 ). 


Proof. Note that u and u 7 satisfy equations (2.6a) and (2.10a), respectively. Thus, by subtracting 
( 2 . 10 a) from ( 2 . 6 a), we obtain that 

f |Vu| p (Vu, Vw) dx— f |Vu 7 | p (Vu 7 , Vu) dx+ f (w ,X7v)dx— f (w 7 , Vv) dx = 0, Vu £ Wg’ p (fi), 
Jq Jq Jq Jq 

which, by choosing v :=u — u 7 , yields that 


f (|Vu| p Vu — |Vu 7 | p Vuy, V(u — u 7 )) dx = f (w 7 — w, V(u — u 7 )) dx, Vu € Wg’ p (fi). (2.11) 
Jn Jn 


Next, by following [12] Th. 2.5], we establish the following pointwise bounds for (w 7 —w, V(u — u 7 )) 
in the four disjoint sets: A PI A,, „4, H X 7 , Ay fl T and X 7 fl X. 


.4n A, 

^inXv 


i 7 ni 

i 7 ni 


((w 7 — w)(x), V(u — u 7 )(x)) < 0. 
((w 7 - w)(x),V(«-it 7 )(s)) < 7 _1 0 2 - 
((w 7 — w)(x), V(u — u 7 )(x)) < 0. 
((w 7 - w)(x),V(u-u 7 )(x)) < r y~ 1 g 2 . 


( 2 . 12 ) 


Since, A n Ay, A PI X 7 , Ay PI X and X 7 n X provide a disjoint partitioning of fi, (2.11) and the 
estimates in ( 2 . 12 ) imply that 


(|Vu| p Vu — |Vu 7 | p Vii 7 , V(u — u 7 )) dx < / 7 g 2 dx. 


(2.13) 


Next, we divide the proof in two cases: p > 2 and 1 < p < 2. 











p > 2: In this case, 138, Lem. 2.1] implies the existence of a positive constant C p , depending 
on p, such that 

(|Vu(.x)| p V : u(x) — |Vu 7 (x)| p Vu 7 (x), V(u — u 7 )(x)) > C p \Vu(x) — Vu 7 (x)| p , a.e. in 1). 


Therefore, by plugging the inequality above in (2.13), we have that 


C p / |Vu — Vu 7 | p dx < / 7 1 g 2 dx, 

J si Jq 


which implies that 


1 ^ II w^ ,p ^ 


o 


g 2 meas(n) \ p 

C P l 


(2.14) 


Finally, since 12 is bounded, (2.14) allows us to conclude that it 7 —> u strongly in ll / n 1 ’ p (0), for 

p > 2. 

1 < p < 2: In this case ( 351 Lem. 2.1] implies the existence of a positive constant D p , depending 
on p, such that 

(| Vu(x)\ p 'Vu(x) — | Vu 7 (x)| p Vu 7 (x), V(u — u 7 )(x)) > D p 1^^L , a.e. in 12. 


(|Vu(x)| + |Vu 7 (x)|) 2 -p’ 


Thus, if we consider this inequality in ( |2. 13 ), we have that 

r |Vu — Vu. |2 


Dn 


oi 


— dx < / 7 1 g 2 dx. 


P 7o(|Vu| + |Vu 7 |)2-p , n 
On the other hand, note that Holder’s inequality implies that 


/ Q |V(u-u 7 )| p dx = f n - 


p(2—p) 

(|VlZ| + |Vu 7 |) 2 


-.(|Vu| + |Vu 7 |)^dx 


< 




(|Vu| + |Vu 7 |) 2 -p 


[/n(|V«| + |Vu 7 |) p cfa] 


2-p 


This last inequality and (2.15) yield that 


.._ , , p (q 2 meas(Q)\ p 

|V(u — w 7 )| p dx ) — y — 


(|Vu| + |Vu 7 |) p dx 


2 -P 

2p 


which implies that 


(2.15) 


Oil W, 


0 


g 2 meas(Q)\ p 
d p1 


ol W 1 ’ 1 " "L 11 "O' II w 1 


'711 W p ’ p 


2-p 

2 P 


Finally, since the sequence {u 7 } is bounded in Wq 1,p (12) (see Lemma 2.5), there exists a positive 

D p meas(Q) 


constant D p such that 


\ u U 7llw [ }’ p — 


1 1 
ryp 


(2.16) 


which, since 12 is bounded, allows us to conclude that it 7 —> u strongly in W 0 1 ,p ( 12 ). □ □ 


9 

























3 Preconditioned Descent Algorithms 


In this section we analyze the application of descent algorithms for solving the regularized problem 


(2.8). We divide this study in two cases: 1 < p < 2 and p > 2. Thus, we need to consider all the 


particular issues that arise in these two scenarios, such as existence of admissible descent directions. 

Descent methods work by finding, at the current iterate u k G V, a search direction w k € V such 
that J-y(uk + tw k ) is decreasing at t = 0, i.e., such that 


(Jy(uk), Wk)v*,v < 0 . 


Here, V stands for a Banach space and V* for its dual space. Although this kind of algorithms are 
usually suitable for differentiable problems, several issues arise. Mainly, the descent provoked in the 
function can be very small. This problem usually appears when the contour maps of the functional 
are very prolonged near the minimizer. Further, in the particular case of problem (2.8), since this 
problem involve the p-Laplacian operator, the difficulties associated to this structure need to be 
taken into account (see ESI). 

An innovative idea to deal with these issues is to use a suitable preconditioner in the computation 
of the search direction. In [25] . the authors successfully implement this idea in a finite dimension 
setting for the p-Laplacian problem. Here, we propose and analyze a similar approach in function 
spaces, for the regularized problem (2.8). In fact, we determine the search direction w k by solving 
the following equation 

Pk{wk,v) = -{J'Ju k ), v) v *,v, Vu G V, 


where V is a suitable Banach space and the form P k : V X V —> M is chosen as a variational 
approximation of the p-Laplacian operator. 

By taking into account the last discussion, we obtain the following general algorithm. 


Algorithm 3.1. Initialize uq G V and set k = 0. 
For k = 1, 2,... do 


1. If ,7( (u fc ) = 0, STOP. 

2. Solve P k (w k ,v) = —(J^(u k ), v)v*,v, Vu G V, for a descent direction w k . 

3. Perform a line search algorithm to determine the step size a k . 

4. Update u k+ \ := u k + a k w k and set k = k + 1. 


Several issues arise when discussing the convergence properties of this algorithm. By considering 
the discussion in |24(, Sec. 2.2.1], global convergence of this algorithm depends on the admissibility 
of w k and a k . 

Admissibility of search directions w k depends on the way in which we define P k . Thus, since 
the behaviour of J 7 depends on the value of p, existence and admissibility of descent directions will 
be discuss in the next sections considering the two cases 1 < p < 2 and p > 2, separately. 

On the other hand, the line search strategy in step 3 of Algorithm |3.1| can be performed in 
several ways. Exact line search algorithms, i.e., algorithms which find a k such that 


J 7 (u fc + a k w k ) = min J 1 {u k + aw k ), 

O !>0 


are known to be expensive, specially when the iterate is far from the solution [40] . Therefore, we 
will use inexact line search techniques. Further, in order to proof convergence for descent algorithms 
like 3.1 these inexact techniques need to be efficient, according to the following definition. 
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Definition 3.2. A line search strategy is called efficient if there exists a constant £ > 0, independent 
of Uk and Wk, such that 


J'yi'U'k d" Oik) — C 


(jy(Uk) } W k)w~l’P' ,Wg’ P 




2 


A classical line search strategy is the so called Wolfe-Powell rule, 
accepting a positive steplength a k if 


This method consists in 


J 7 (u k + aw k ) < J'yfuk) + oiak(J!y(uk), w k )v*,v, (3.1a) 

( J!y(u k + aw k ) , w k ) v *,v > o- 2 (J^(nfc) , w k )v*,v , (3.1b) 

where 0 < o± < 02 < 1. Wolfe-Powell rule is known to satisfy the previous efficiency requirements 
and it will be used as a central requirement in the coming convergence results. 


3.1 The 1 < p < 2 case 


In this section, we construct an algorithm, based on Algorithm 3.1, for the problem (2.8), when 
1 < p < 2. Due to the structure of the problem, we will analyze this case in function spaces. 
Therefore, we discuss the space V in which the algorithm is constructed, define the bilinear form 
P k (-, ■), analyze the equation Pk(w,v ) = — (J 7 (u), v)v*,v, and, finally, we write the algorithm and 
prove a global convergence result. 

Definition 3.3. Let l<p<2, e>0 and u € Wo’ p (D). We define Hq (fi) as the completion V(D) 
with respect to the norm 


I Z \\HS = 


(e+ \Vu\) p - 2 \S7z\ 2 dx 


Theorem 3.4. Let 1 < p < 2, e>0 and u E Wq’ p (D). Then, Hq(Q) is a Hilbert space with the 
inner product 

{z,w) H % = [ (e + \Vu\) p - 2 (Vz,Vw)dx. (3.2) 


Furthermore, the following inclusion holds, with continuous injections 

H^D) C H%(D) C W 0 1 - P (fi). 


(3.3) 


Proof. Let us start by pointing out that (3.2) is a positive definite bilinear form, which fits the 
structure analyzed in [9j p. 214] and [121 pp. 268-269]. 

Next, we analyze the coefficient (e + |V«|) p-2 . First, note that 


(e + \Vu(x)\y~ 2 = 


(e + |V'u(.t)|) 2_ p 


, a.e. in D, 


which implies, since 2 — p > 0, that 


< 


a.e. in Ll. 

) 7 


(e + |Vw(x)|) 2- p e 2- P’ 

The last two expressions yield that 

(e + |Vu|) p - 2 EL°°(D) Cl 1 ^). 


(3.4) 
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Now, note that 


[(e + |Vu(x)|) p 2 ] 1 = (e + |Vu(x)|) 2 p , a.e. in fi. 


Since u £ Wp ,p (fi) and 2 — p < p, we can state that 


{e + \Vu\) 2 ~ p £ L\n). 


(3.5) 


Consequently, (3.4), (3.5), |9j Lem. 3.3] and [321 P- 268-269], yield that the Hilbert space iLg(H) 
is well defined. 

We now prove (3.3). Let z £ Hq(Q). First, note that, thanks to (3.4), there exists a positive 
constant C\ > 0, such that 


(e + \Vu\) p ~ 2 \Vz\ 2 dx < Ci [ \Vz\ 2 dx, 

i Jn 


which implies the existence of a positive constant C\ such that 

\\ z \\h% < Ci\\z\\ H i. 

Further, let z £ Hq (f2). Holder’s inequality implies that 

|VH P . . P(2~P) 

-^(e + |Vu|) 2 dx 


(3.6) 


Jn\Vz\ p dx = f a - p(2 _ p) 

(e + |Vu|) 2 


< 


ft 


|V*| : 


n (e+ |Vu|) 2 -p 


dx 


[J n (e + \Vu\) p dx]^ P 


Next, since u £ H / 0 1 ’ P (H), the last expression implies the existence of a positive constant C 2 such 
that 


\Vz\ p dx < C 2 


(e+ \Vu\) p ~ 2 \Vz\ 2 dx 


which implies the existence of a positive constant C 2 such that 


z \\w}* - C ’ 2 W Z \ 


HU 


(3.7) 


□ 


Summarizing, (3.6) and (3.7) imply that 

77q(H) C H%(n) C W 0 1 ’ p (fi), 

with continuous injections. □ 

We propose our algorithm, considering that V := Hq(CI), for some suitable u £ H , ’ 0 l ’ p (O). 
Moreover, it looks natural that the form P' k will be defined as follows. 

P k (w,v) := [ (e + \Vu\) p ~ 2 (Vw,Vv)dx. 

Jn 

Here the small parameter e > 0 helps the algorithm to handle possible degeneracy when Vu = 0. 
Note that P k is a linearization of the weak form f n | \7u\ p ~ 2 (Vu, Vu). 

Next, note that Hq (H) C Wg 1,p (r2), for all u £ H / 0 1,P (H). Next, we define by J 7 (u) the restriction 
of J!y(u) to Hq(Q). Therefore, we can state that J' 7 {u) £ Hq(Q)* and that 


(eXy(u), v) jju* uu (J 7 (u) , r)jy_ 1( p/ Vv £ Hq (H). 


(3.8) 
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For further details, we refer the reader to |41 Rem. 3, p. 136]. 

Summarizing, we need to analyze the following variational equation 

[ (e + \Vu\) p - 2 iVw,Vv)dx = -(J'(u),v) H ^ H u, Vr £ (3.9) 

It is clear that a solution for a similar equation will play the role of the descent direction in our 
Algorithm. Therefore, we need to prove that this equation has, at least, one solution in Hq (fi). 

This existence result is a direct consequence of the Riesz-Frechet representation theorem (see 
[U Th. 5.5]). In fact, we know that Hq(Q) C Wq' p (Q.) is a Hilbert space for 1 < p < 2. Moreover, 
we know that Jq(c + |Vn|) p ~ 2 (Vu;, Vu) dx is the scalar product of this Hilbert space. Consequently, 
the Riesz-Frechet representation theorem implies the existence of a unique w £ Hq(Q) such that 


[ (e + | Vri|) p_2 (Vtc, Vv) dx = -{J'Ju ) , v) H u* H u, Vu £ Hq(Q). 

Summarizing, the Algorithm |3.1| takes the following form for 1 < p < 2. 

Algorithm 3.5. Initialize uq £ W^ ,P (Q) and set k = 0. 

For k = 1, 2,... do 

1. If J'^Uk) = 0, STOP. 

2. Find a descent direction w k £ i?Q fe (f2) by solving the following variational equation 

f n (e + | Vu k \) p ~ 2 (Vwk, Vv) dx = -(J'(u k ), v) H " k * H u k 

° ’ ° (3.10) 

= ~ In \Vu k \P~ 2 (Vu k ,Vv) dx - g 7 f n J^^i) dx + / n / v dx , \/v £ i?£ fc (fi). 


3. Perform an efficient line search technique to obtain a k - 

4. Update u k +i '■= u k + a k w k £ Wo’ p (fi) and set k = k + 1. 

Clearly , the equation (3.10) has a unique solution Wk £ HQ k (Q) C IUq'^O), for all k £ N. Thus, 


Algorithm 3.5 is well defined. However, it is mandatory to prove that w k £ ITq ’ p (Q) is, indeed, an 


admissible descent direction. First, we prove that Wk is a descent direction. In fact, note that from 


(3.8) and (3.10), we obtain that 


(Jryiu) , Wk) w -i, P ' t wi ,p — (Jiy (u), Wk)i 


= f n {e + \Vu k \) p 2 \Vw k \ 2 dx 
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which yields that 

(J^(u ), jyUp < 0. (3-11) 

Next, let us discuss the admissibility of w k - Note that if we had defined Pk as the variational 
version of the Laplacian operator, i.e., P k (u,v) = f n (Vu, Vv)dx, the sequence generated by the 
associated version of the Algorithm 3.5 would be such that {u k } C C H / 0 1,P (H). In this 


case, it is possible to state the existence of q < 2 such that Hq(Q.) C IUq’^O) C Wo’ p (fi) (see 
HSU)- On the other hand, note that the sequence {u(} generated by Algorithm |3.5| yields that 
{uz} £ U (p) C lF ( j ,p (Q). These arguments suggest that, using interpolation theory HU, a 


similar inclusion result can be obtain for Pk(w,v) := f n (e + |Vu|) p 2 (Vw, Vv) dx. Thus, we make 
the following assumption. 
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Assumption 3.6. There exists q, 1 < p < q < 2, such that U j£N H^(Q) C W 0 ll? (fi) C <*(0). 


Proposition 3.7. Let {uk} be the sequence generated by Algorithm 3.5 and suppose that the step 
length ak satisfies the Wolfe-Powell conditions (3.1). Furthermore, let us suppose that the Assump- 
tion\3.6\ holds. Then, the Zoutendijk condition is verified, i.e., 


OO 

COS 2 4>k = OO, 

k =0 


(3.12) 


where cos 4>k = ~ 


(Jfi u k) ,' w k ) w _ 1 y w l,p 


0 




Proof. First, note that Theorem 2.1 implies that the functional J 7 is bounded below in VFq’ p (Q). 
Next, let us recall that the functional J 7 can be written as 

J 7 (tt) = F(u) + 


where F(u) and C/ 7 (Vii) are given in Section 2.2 It was previously stated that both T and Q 7 are 
continuously differentiable in \V^' :P (ii). Moreover, it is known that T is actually twice differentiable, 
since this functional represents the variational version of the Dirichlet problem for the p- Laplacian 
operator (see mm ). Thus, T has a Lipschitz continuous gradient in W 1,p/ (f l). On the other 
hand, in Section [2.2| we stated that 


<s;(vu) 


v )w~ 1 ’P l ,wF r ~ 9 T 


r (v«, Vv) 

'n max( 5 , 7 |Vu|) 


dx. 


Next, thanks to the Assumption |3.6[ the max function involved in the last expression is slantly 
differentiable (see |2TJ). Consequently, we can state that Q'{S7u) is slantly differentiable in Wq P ({T). 
Therefore, thanks to m Th. 2.6, pp. 1205], Q!y is Lipschitz continuous in W 1,p, (fi). 

Summarizing, we know that J 7 is bounded below and continuously differentiable in 
and its gradient is Lipschitz continuous in W~ 1,p (fi). Therefore, since we assume that ak satisfies 
the Wolfe-Powell conditions, all the hypothesis of Zoutendijk theorem are satisfied (see, for instance, 
m pp. 29] and pLUl Lem. 2.5.6] and the references therein ), and, consequently (]3. 12) holds. □ □ 


Theorem 3.8. Let {life} be the sequence generated by Algorithm \3.5\ and suppose that the step length 
ak satisfies the Wolfe-Powell conditions (3.1). Furthermore, let us suppose that the Assumption 


3.6 holds. Then, the sequence {ttfc} converges to the uniquely determined global minimum of J 7 


Proof. First, note that the Hanner’s inequality m and the convexity of function ip imply that J 7 is 
a uniformly convex functional in Wq’ p (Q). Further, Proposition 


3.7 


guarantees that the Zoutendijk 
condition holds. Therefore, the result directly follows from |18l Th. 4.7]. □ □ 


3.2 The p > 2 case 


In this section, we construct an algorithm, based on Algorithm 3.1, for a discrete approximation 


of the problem (2.8), when p > 2. Our first aim was to construct an algorithm in function spaces. 


However, the structure of the problem prevents us from this goal. Particularly, there are regularity 
issues regarding the search direction. Indeed, we have the following result. 
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(3.13) 


Theorem 3.9. Let p > 2 and <p G W 1,p, (Jl). 


Then, the variational equation 


/ (Vw,X7v)dx = (v?, v) w _ iy w i, P , \/v G W 0 1,p (n) 
has a unique solution w G Rjj ,p (H). Furthermore, there exists K > 0 such that 

K \\M\wl* - Ill’ll w-^r' ^ W w Ww£’ p - ( 3 - 14 ) 

Proof. Since i) C R 2 is assumed to be a bounded domain with regular boundary, [371 Th. 4.6] 
immediately implies the result. □ □ 

Note that J 7 : W 0 1,P (J4) —> M, which implies that J!y(u) G W~ 1,p ' (Lt), for all u G Wg’ p (fi). 
Therefore, it is possible to find a unique solution w k for the following equation 


[ {Vw k ,Vv)dx = ~(J'(u k ), v) w _ hp , w i.p, Vv G Wl' p {yi). 
Jn ’ 0 


However, w k G Wq’ p (fl) D V 
algorithm like Algorithm 
extended to more genera’ 


3.1 


/q ,p (H) for p > 2. This fact prevents us from directly constructing an 
since u k + 1 = u k + afcuifc G W 0 1,p (Q). Moreover, Theorem 


el 


3.9 


can be 


iptic forms than the Laplacian. These results can be found in, e.g., 
Consequently, the regularity issue prevails, for several elliptic choices for P k . 

A possible solution for this issue is to pose the problem in a suitable H s (Ll) space, with sGR 
such that H S (Q) C W 0 1,p (f2). Indeed, it is known that for p > 2 and u G Wo’ p (fi), the following 
inclusions hold, with continuous injections (see 0 ) 


Furthermore, it is possible to state that (see [10, Prop. 1 pp. 96]) 

H S (Q) C H\n), Vs > 1. 


(3.15) 


(3.16) 


Thus, pUl Rem. 2 pp. 96], (3.15) and (3.16) yield 


the existence of a s G M such that 


H S (Q.) C W 0 1 ,p (fi) C H\n). 


Therefore, we can define P k as the scalar product in H s (Ll). However, several technical challenges 
arise with this idea. For instance, the actual value of s is unknown, and the numerical realisation 
of the search direction requires the implementation of the Fourier transform of several functions. 
We consider that all of these issues are beyond the scope of this paper, and will be considered in a 
future contribution. 

Another possible idea to overcome the regularity problem is given by a smoothing step. 

W^ p \n) 3w k ^ W k G Wo llP (fi). 


In Sec. 6], the author discusses the definition and properties of such a procedure. Though 
this smoothing procedures are designed for fixing regularity issues in function spaces like the one 
we have is this paper, they need several technical assumptions. These assumptions, at least in this 
context, can be very restricitve and can even reduce the admissible set of solutions for equation 
P k (iu, v ) = , v ) to the empty set. On the other hand, it is known that in finite dimensional 
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spaces no smoothing step is needed, so we can define Wk as the search direction for the descent 
algorithm (see [l3j Sec. 6.1]). 

By taking into account the argumentation above, we consider that the best solution is to 
analyze the problem with a “discretize then optimize” approach. Thus, we propose a finite element 


discretization of the problem (2.8). Next, we propose and study a preconditioned algorithm for the 


case p > 2 in finite dimension spaces. 

We propose a discretization with first order finite elements, following ideas in [2) 19]. Thus, 
let T h be a regular triangulation, in the sense of Ciarlet, of 0. Next, let Q h be a polygonal 
approximation to hi, given by fl h = UreT h O where all the open disjoint regular triangles t have 
maximum diameter bounded by h. Further, for any two triangles, their closures are either disjoint 
or have a common vertex or a common side. Finally, let {Pj}j=i,...,N be the vertices associated 
with the triangulation T h . Hereafter, we assume that Pj E dQ h implies that Pj E dtt and that 
Q h C H. In this paper we will only consider first order approximation, because of the limited higher 
order regularity for the solutions of the p-Laplacian (see [25] and the references therein). Taking 
the above discussion into account, we introduce the following finite-dimensional spaces associated 
with the triangulation T h 

Wq := {u E C(VL h ) : v\ T G Pi,Vr G T h and v = 0 on dQ h }, 

where Pi is the space of polynomials with degree less than or equal to 1. 

Thanks to these defintions, we can introduce the following finite element version of the problem 


( 2 . 8 ): 


min J^(u h ) := - 
u h £Wfr 


:= — f \Vu h \ p dx + [ u h ) dx — [ f u h dx. (3-17) 

V J ci h Jn h J n h 


Theorem 3.10. Problem (3.17) has a unique solution u h G Wq. 


Proof. This result is a direct consequence of the fact that Wjf is a closed subspace of Wq’ p (Q) (see 
PH Sec. 3.2]). □ □ 

As stated in the previous section, in finite dimensional spaces is not mandatory to use smoothing 


steps to construct preconditioned descent algorithms for problems like (3.17). Further, in this case 
we know that (see mm) 

C W 0 1,P (Q) C (3.18) 

Thanks to this fact, we can consider Wq a Hilbert space with the norm induced by 77 q(H), which 
we will note by || • \\ W h. 


Summarizing, we propose the following algorithm for problem (3.17) with p > 2. 
Algorithm 3.11. Initialize uq G Wq and set k = 0. For k = 1,2,... do 

1. If 4'(?4) = 0, STOP. 

2. Find a search direction by solving the following variational equation 

fn( Vw k’ Vv ) dx = -( J y( u k) > v )(w£)*,w£ 

= - Sa t 2 ( VujS, Vv) dx - 97 /n dx + J a J v dx, V„ € Wf. 

3. Perform an efficient line search technique to obtain a*,. 

4. Update u£ +1 := + ak w k an d set k = k + 1. 


(3.19) 
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Proposition 3.12. The equation (3.19) has a unique solution € Wq . Furthermore, this solution 
w l f is an admissible descent direction for J 7 (uk), i.e., it satisfies that 


{Jy ( u k ) ) W k)w!t* Wn ^ ^ ^) 


and the following admissibility condition 

{Jy ( u k ) j W k)wf*,w/t 


\ w k\\w$ 


0 ^ 0 => || J 7 (^fc)llwh* > 0. 

«—>■ OO u K—>■ OO 


(3.20) 


Proof. Existence of a unique solution directly follows from the fact that Wjf is a Hilbert subspace 
of Hq(SV) with the induced norm of this space. Therefore, wi is the Riesz representation of the 
functional —J 7 (wj£) in the space Wq (see j25l Sec. 3.1]). Furthermore, thanks to (3.18), from 
(3.19) we can conclude that 


(Jy ( u k ) > W k)wh*,wh ~ ~\\ w k\\'wf < ^ e 


h II2 


(3.21) 


which yields that w% is, indeed, a descent direction for jfj' (u(i). Finally, since is the Riesz 


representation of J^' (u%) in Wq, we have that 


T ~ II Jy ( u k)\\ W f*- 


This last identity, together with (3.21), yield that 

jh' f n .h 


which immediately implies (3.20) 


(Jy ( u k) > W k)wft*,w£ — W^y ( U k ) 11 Wfr* 11 W k 11 W$ ’ 

□ 

Then, 


Theorem 3.13. Let w£, otk and ui generated by Algorithm 


3.11 


lim J 7 (u^) = 0. 


k—¥ oo 


□ 


(3.22) 


Proof. Since vfi: satisfies (3.20) and aik is calculated by an efficient line search algorithm, admissi¬ 
bility of these two sequences is guaranteed. Therefore, since all the hypothesis of [23 , Th. 2.2] are 
fulfilled, we can conclude the proof. □ □ 


4 Numerical Implementation 

In this section we discuss all the issues related to the numerical implementation of the algorithms 
developed in the last section. Further, we present several numerical experiments to show the 
behavior of these algorithms. Such experiments are concerned with the two cases analyzed during 
this paper: 1 < p < 2 and p > 2. The case p = 2 has been widely analyzed, by using a similar 
regularization approach, in pa na nu. Moreover, we focus our experiments on the numerical 
simulation of the laminar flow of a Herschel-Bulkley fluid in a pipe. Therefore all the experiments 
have been carried out for a constant function /, which represents the linear decay of pressure in 
the pipe. 
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4.1 Discretization issues 


In this section we describe the finite element implementation that we use in all the numerical 
experiments. Let us start by pointing out that we use the same finite element approach described 
in Section 3.2 Thus, we recall the finite dimension space 


:= {v E C(Q h ) 


E Pi,Vt E T h and v = 0 on dtt h }, 


where Pi is the space of polynomials with degree less than or equal to 1. We note the basis functions 
of Wq by <pj, j = 1,... ,n and we assume that card(T h ) = m. Further, we use the notation it for 
the coefficients of the approximated functions u h . 

By following ideas in [121 Sec. 4], we use the following discrete version of the gradient 


X7 h : = 



E M 


2 mxn 


(4.1) 


where 8* := ^ and 8% := 

T k 


d<Pi(x) 


dx 2 


are the constant values of 


Tk 


8x2 

8<fi(x) 

dx\ 


for i = 1 ,...,n and r k E T h . Note that 


dipi(x) 


Tk 


dxi 


and 


Tk 


and in each triangle r k , respectively. Consequently, 


\7 h lt is the approximation of \7u h (x). 
Next, let us introduce the function £ : 


given by 


£(w)k = \{vJ k ,w k+m )\ T , k = l,...,m. 

Therefore, we calculate \Vu h (x)\ by £( X7 h lt ). Note that represents the value of | SJu h (x)\ 

at each triangle r k E T h . 

Finally, we discuss the implementation of J fi (e+ |Vu|) p_2 (Vui, Vv) dx, | Vu\ p ~ 2 (Vu, Vv) dx 

and 57 |) dx ‘ using the Galerkin’s method, we obtain the following 

• fn( e + \Vu\)P~ 2 (Vw,V(pj)dx « E"=i w i ErkGT^ f Tk (e + dx, 


• /ol Vu l P 2 (Vu, X7(pj)dx « YJi=i u il2rk&T h / T J£(V ft -h > ) fc ) p 2 (V^,V^)dx and 


fn 91 


(V'U.V^j) , 


. ■sr^n 




for j = l,...,n. Next, note that the terms (e+^(V h lt)k) p ~ 2 , (£,(V h lt)k) p - 2 and max(< 7 , 7 £(V h l/ , )fc) 
are constant at every triangle r k . 

Thus, by using ideas in [[5], we obtain a matrix approximation A^ u E M nxn , E M nxn and 

Ajmax £ M nxn , for any of the forms in the expression above. The entries of these matrices are 
given by 


• (a e ,u)i,j = Er fc eT'‘( e + C(V^)A : ) p 2 f Tk (V< Pi ,V< Pj )dx, 


• Mzj = 9lY,T k &T h ^^y h ~^)k) p 2 4(V^,V^)d*and 


• (a u ,max)ij - 57 Yj Tk &T h ma x(g,rt(V h lt) k ) V ^i) dx ~ 

Finally, by following ideas in [2], we approximate the right hand side as follows 



,n, 
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where the quadrature rule Q T is given by 


1 3 

Q r (v) = -meas(r) v(ai), with ai, i = 1,..., 3 the vertices of r S T h . 

3 i =l 

Remark 4.1. It is remarkable to point out that due to the proposed structure, the Algorithms 


3.5 and 3.11 only need to solve one linear system at each iteration. In fact, Algorithm |3.5| and 


Algorithm 3.11| require the solution of linear systems like 


A h e , u ^ = V f, and A h ^ = r,!j, 

respectively. Here, A^ u is given above, A h is the classical stiffness matrix and and rjfy are the 


F.E.M. approximation of the right hand side of equations (3.10) and (3.19), respectively. Note that 


matrix A^ u depends on Uk, but does not depend on w k . Further, A h does not depend neither on Uk 
nor in w k . This fact implies that the linear systems can be easily solved by any direct or iterative 
method and does not represent a large computational effort. 


Remark 4.2. (Stopping Criteria) We stop the Algorithms 3.5 and 3.11 as soon as the expression 
| is reduced by a factor of 10 -6 . Here ./(, h (~etk) stands for the FEM discrete version of 
J'(u k ) and is given by 

:= k + ^u,max~^fc — 7- 

This stopping criteria are popular for steepest descent algorithms, since they are easy to implement 
and provide enough information about the convergence behavior of the algorithm (see [293). 


4.2 Line search algorithms 

As stated in Section [3| we need to focus on the implementation of efficient inexact line search 
methods. One typical technique is the backtracking line-search algorithm. The general idea behind 
this approach is to take ak = 1. Then, if Uk + a k w k is not acceptable, in the sense that a descent 
condition on J 7 is not fulfilled, a k is reduced (“backtracked”) until u k + a k w k is acceptable. 

We propose to use an algorithm which uses polynomial models of the objective functional for 
backtracking, which is detailed in PH Sec. 6.3.2], In this section, we briefly describe this algorithm 
and the main ideas behind it. 

The central discussion in a backtracking algorithm is how to reduce a k . Usually, the back¬ 
tracking algorithm is implemented by taking ak = -^a k , so ak is reduced to half at each iteration. 
This procedure can be inefficient since usually needs several iterations to achieve convergence, and, 
moreover, the step sizes can be very small. 

In this paper, following ideas in [TB1 Sec. 6.3.2], we propose a reduction strategy for a k based 
on polynomial models of the objective function. 

Let us start by introducing the following function 

<y9 fc (a) := J 7 (u fc + aw k ). 

Next, by using the current information of J 7 , we take a k as the approximation to the value that 
minimizes tp k (a), Le., a k ~ argmin ip k (a). 

First, note that the following information about <p k is available. 

Wfc(0) = J{u k ) and ^(0) = (J 7 (u fc ), w k ). (4.2) 
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Further, once we calculate J^(u k + w k ), we know that 

^fc(l) = J-y(Uk + Wk). (4.3) 

Next, if Jj(u k + Wk) does not satisfy the descent condition (he., </?&(!) > ^(0) + a\ip' k (ti)), we 


construct the following quadratic model for tpk by using (4.2) and 


m 2 {a) := (<p k ( 1) - <p k (0) - <ff k (0))a 2 + <p' k {0)a + tp k { 0). 


It is easy to prove that 


a 2 = 


-^4(°) 


(4.4) 


2 (^( 1 ) -^( 0 )-^( 0 )) 
is a stationary point of m 2 , i.e., it satisfies that m' 2 (a) = 0. Moreover, we have that 

m 2 (a) = 2(<p k (l) - ip k { 0 ) - ^ fc ( 0 )) > 0 , 

since ^(1) > p k (0) + <Ji<p' k (Q) > ip k { 0) + <// fc (0). Thus, we conclude that a 2 minimizes the model 
m 2 and, since ip k ( 0) < 0, we have that a 2 > 0. Consequently, we take a k := a 2 . 

Now, since <^fe(l) > <y?fc(0) + oap' k { 0), from (4.4), we have that 


a 2 < 


2(1 — ci) 


Therefore, (4.4) gives an implicit upper 


This fact implies, provided ip k { 1) > ip k { 0), that a < 
bound of ~ 7 for a 2 on the first backtrack. On the other hand, if <PkiX) >> <^(0), a can be very 
small. This fact suggests that <pk{oi) is probably poorly modeled by a quadratic function is this 
region. In order to avoid too small steps, we impose a lower bound of yy. Therefore, if at the first 
backtrack at each iteration we have that cb < 0 . 1 , the algorithm next tries a k = yy 


Now, suppose that <^( 0 : 2 ) does not satisfy (3.1a), which implies that we need to backtrack again. 
In this case, we have the following information available: ip k {0) = J{uk ), ^4(0) = (J'(u k ), w k ) 
and the last two values of p(a). Therefore, we use a cubic model of ip fitting all these pieces of 
information, and set a k to be the minimizer of this new model. This procedure is justified since a 
cubic polynomial can perform better when modelling situations where J 7 has negative curvature, 
which are likely when (3.1a) is not achieved for two possible values of a (see m pp - 128]). 

The construction of this cubic model is as follows. Let a p and a 2p be the last two previous 
values of a k - Then, the cubic that fits ip k ( 0), <p' k { 0), ip k {oi p ) and p k (& 2 p) is given by 


m3 (a) := ca 3 + da 2 + </4(0)a + ^fc(O), 


where 


1 


Otp &2p 


KA-p 

-&2p 

Op 

p 


( <fk{a P ) ~ Pk{ o) - ^fc(0)a p \ 
V^fe(«2p) - <£fc(0) - ip k (0)a 2p .J 


Further, it is easy to prove that the minimizer of m 3 is 


«3 = 


—d + y / d 2 - 3cv?' fc (0) 


3c 


(4.5) 


In (16] is established that if <p(a p ) > 99 ( 0 ), then 0:3 < |a p , but this reduction is considered too 
small. Therefore, we impose the upper bound b = 0.5, which implies that if 0:3 > \ol v , we set 
a k = \ol p . Also, since 0:3 can be an arbitrarily small fraction of a p , we again impose the lower 
bound a = yy , i.e., if 0:3 < yy a p , we set a k = yy a p . 

Summarizing, we have the following line search algorithm. 
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Algorithm 4.3. Let 07 e (0, \) and set «o = 1. 

1. Decide wheter J 1 {uk + otk) > Jj{uk) + a\ak{Jly{uk ), Wk) holds. If so, STOP and set otk = ao- 
If not: 


2. Decide wheter steplength is too small. If so, STOP and terminate algorithm: routine failed 
to locate satisfactory Xk +1 sufficiently distinct from Xk- If not: 

3. Decrease cc by a factor between 0.1 and 0.5 as follows: 

(a) On the first backtrack: set a*. := a 2 = arg min m 2 (a:), but constrain the new to be 

> 0.1. 

(b) On all the subsequent backtracks: set ak := 0:3 = arg min m 3 (cc), but constraint the new 
a*, to be in [0.1 a p , 0.5a p ]. 

4. Return to step 1. 


Here, the parameter 07 is set quite small, usually in the order of 10 4 . Further, (4.5) is never 
imaginary if 07 is less than t (see [16| pp. 129]) 


Note that this algorithm only implements the first Wolfe-Powell condition (3.1a). The curvature 
condition (3.1b) is not usually implemented because the backtracking technique avoids excessively 
small steps. It is established that the bounds in the algorithm on the amount of each calculation 
of a make the curvature condition to hold (for further details and examples see [I~6l pp. 126-129] 
and the references therein). 


4.3 Numerical Results: Case 1 < p < 2 


In this section, we focus on the behavior of Algorithm |3.5| In the next experiments, we consider 
that the problem (2.8) represents the flow of a Herschel-Bulkley fluid with 1 < p < 2, so we are 


in the case of a shear-thinning material. Further, we consider a constant /, which represents the 
linear decay of pressure in the pipe. In this context, the constant g plays the role of the plasticity 
threshold and it is modelled by the Oldroy number (see [25]). For further details in the mechanics 
of these problems, we refer the reader to m na ESI and the references therein. 

Hereafter, we use uniform triangulations described by h, the radius of the inscribed circumfer¬ 
ences of the triangles in the mesh. In the next examples, we use the values 7 = 10 3 and e = IIP 6 , 
and we initialize the algorithm 3.5 with the solution of the Poisson problem — Au!q = f h . Further, 


we stop the algorithm by using the stopping criteria described in Remark |4.2 


4.3.1 Experiment 1 

In this experiment, we set fl C M 2 to be the unit ball, and we compute the flow of a Herschel- 
Bulkley material with p = 1.75. We analyze the behavior of the algorithm with g = 0.2 and / = 1, 
and we use a mesh given by h ~ 0.0086. 

The resulting velocity function and the velocity profile along the diameter of the pipe are 
displayed in Figure [2] The graphics illustrate the expected mechanical properties of the material, 
i.e., since the shear stress transmitted by a fluid layer decreases toward the center of the pipe, the 
Herschel-Bulkley fluid moves like a solid in that sector. This effect explains the flattening of the 
velocity in the center of the pipe. 

In Table [lj we show the number of iterations that Algorithm |3.5| needs to achieve convergence. 
We also show the value of the value of Jy : hC^k)i the value of the step 07 . and the 
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Figure 2 : Calculated vleocity u (left) and velocity profile along the diameter of the pipe (right). 
Parameters: p = 1.75, g = 0.2, 7 = 10 3 and e = 10" 6 . 


it. 

l^iWI 

k) 

Oik 

l.s. it. 

1 

2.070e-3 

-0.022393 

1.0000 

0 

2 

8.758e-3 

-0.027232 

0.4199 

1 

3 

2.959e-3 

-0.028522 

0.3390 

1 

4 

5.582e-4 

-0.028778 

0.2600 

1 

5 

3.532e-4 

-0.028911 

0.1231 

2 

6 

7.549e-4 

-0.029028 

0.1864 

1 

7 

6.590e-4 

-0.029057 

0.0788 

2 

8 

4.865e-4 

-0.029091 

0.0558 

2 

9 

1.179e-4 

-0.029101 

0.0485 

2 

10 

6.655e-7 

-0.029107 

0.0424 

2 


Table 1: 
e = 10 -6 . 


Convergence behavior for Algorithm 


3.5 


Parameters: p 


1.75, g = 0.2, 7 = 10 3 and 



Figure 3: Calculated residual | h {uk)\/\J!y h (ito)| for: Algorithm |3.5| (left) and Wolfe-Powell (right). 

Parameters: p = 1.75, g = 0.2, 7 = 10 3 and e = 10~ 6 . 
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e 

it. 

J{u) 


le-4 

10 

-0.029107 

1.114e-6 

le-5 

10 

-0.029107 

7.064e-7 

le -6 

10 

-0.029107 

6.655e-7 


Table 2: Dependence on e for Algorithm 3.5 Parameters: p = 1.75, g = 0.2 and 7 = 10 3 


5 = 0.1 

hi 

h 2 

h 3 

Iter. num. 

9 

9 

9 


1.147e-6 

1.667e-6 

1.661e-6 


-0.0395 

-0.0414 

-0.0416 

to 

II 

O 

to 

hi 

h 2 

h 3 

Iter. num. 

9 

8 

8 


1.490e-6 

7.059e-7 

3.976e-6 


-0.0217 

-0.0231 

-0.0233 

CO 

0 

II 

0 ^ 

hi 

h 2 

h 3 

Iter. num. 

18 

19 

19 

\ J '-y,h( u k)\/\J!y, h (. u o)\ 

4.232e-6 

4.393e-6 

1.342e-6 


-0.0105 

-0.0115 

-0.0116 


Table 3: Convergence behavior for Algorithm 3.5 Parameters: p = 1.5 and 7 


10 3 


number of inner iterations needed by Algorithm 4.3 As expected, the value of the functional is 
monotonically reduced at every iteration. The residual behaves typically as in a steepest descent 
algorithm, as shown in Figure[3j However, | h (u k)\/\J!y /i("^o)| decays faster in the last iterations. 
This fact suggests, at least experimentally, that this algorithm has a fast local convergence rate. 
The step ck*, is also monotonically decreasing and the line search Algorithm |4.3| needs no more than 
two inner iterations to calculate the step. 

Finally, in Table [2] we compare the behavior of the Algorithm 3.5 for different values of the 
parameter e. It is clear that the performance of the Algorithm is similar in the three cases shown. 
Some small improvement can be seen, though, for small values of e. 

Let us emphasize that our method requires a low computational effort to produce results which 
are in good agreement with previous contributions (e.g.,[26|). In fact, we only need to solve one 
linear system per iteration and the line search strategy needs two iterations in average. 


4.3.2 Experiment 2 

In this experiment, we set D to be the unit square (0,1) x (0,1), and we compute the flow of a 
Herschel-Bulkley material given by p = 1.5. We fix / = 3, and we focus on the behaviour of the 
algorithm in different meshes, since we are interested in showing, at least numerically, the mesh 
independence of our algorithm. It is known that smaller values of p imply that the functional 
loses regularity, making the problem a bit more challenging. In fact, as g grows, the contribution 
of the less regular component of the functional f Q i/i 7 (Vu) dx increases. This fact complicates the 
numerical approximation of the problem. Therefore, we test our algorithm with several values of g 
to show the versatility of our approach. 

In Table [3j we present the main features of Algorithm 3.5 for several values of g and different 
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Figure 4: Calculated residuals || J / (ufc)||, in different meshes, for p = 1.5 and g = 0.1 (left), g = 0.2 
(center) and g = 0.3 (right). Parameters: 7 = 10 3 . 


mesh sizes: h\ ~ 0.0133, /12 ~ 0.0047 and ~ 0.0029. As expected, the number of iterations that 
the Algorithm needs to achieve convergence increases as g does. However, for a given g , the number 
of iterations is very stable as the mesh size decreases. Also, the evolution of | J' h (ltk)\/\J^ h {u o)| is 
quite similar at every mesh, as shown in Figure [4j These facts show the robustness of our approach 
and numerically verify the mesh independence of the algorithm. 

The resulting velocity functions and the velocity profiles along the diameter of the pipe are 
displayed in Figure[5| As in the previous case, the shear stress transmitted by a fluid layer decreases 
toward the center of the pipe which provokes the solid-like movement in that sector. Further, it is 
expected that if the value of g increases, the flow tends to slow down and the flat zones tend to be 
bigger. This is clearly shown in the figures depicted, which are in good agreement with previous 
contributions (e.g., 


4.4 Numerical Results: Case p > 2 


In this section, we focus on the behavior of Algorithm |3. 11 In the next experiments, we consider 
that the problem (2.8) represents the flow of a Herschel-Bulkley fluid with p > 2, so we are in the 


case of a shear-thickening material (see [26]). Further, we consider a constant /, which represents 
the linear decay of pressure in the pipe. As in the previous section, the constant g plays the role of 
the Oldroy number. For further details in the mechanics of these problems, we refer the reader to 
0Q21I25] and the references therein. 

with the solution of the Poisson problem — Aiig = f h , and we 


We initialize the algorithm 


3.11 


terminate the iterations according to the stopping criteria described in Remark |4.2| 

As in the previous section, we use uniform triangulations described by h, the radius of the 
inscribed circumferences of the triangles in the mesh. 

It is remarkable to state that the classical p-Laplacian problem (ie., (2.8) with g = 0) is difficult 


to solve when p + l/(p — 1) is large. This issue needs to be take into account in our case too (see 

E5])- 


4.4.1 Experiment 1 

In this experiment, we set fl to be the unit square, and we compute the flow of a Herschel-Bulkley 
material with p = 4. We analyze the behavior of the algorithm with g = 0.2 and / = 3. We work 
with a mesh given by h ~ 0.0029, and we use the value 7 = 10 3 . 

The resulting velocity function and the velocity profile along the diagonal of the square pipe are 
displayed in Figure [ 6 ] The graphics illustrate the expected mechanical properties of the material: 
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Figure 5: Calculated u for p = 1.5 and g = 0.1 (top, left), g = 0.2 (top, right) and g = 0.3 (down, 
left). Velocity profile for the calculated velocities along the diagonal of the square (down, right). 
Parameters: 7 = 10 3 . 




Figure 6 : Calculated u for p = 4 (left) and velocity profile along the diagonal of the pipe (right). 
Parameters: g = 0.2 and 7 = 10 3 . 
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Figure 7: Calculated residuals \ J' h {uk)\/\J' ry fe (rto)|, for p = 4 and g = 0.1. Parameters: 7 = 10 3 . 


it. 

K^TOI/Kh^o)! 



l.s. it. 

1 

5.4526e-3 

-0.17950 

1.0000 

0 

2 

2.9779e-4 

-0.18069 

0.4974 

1 

3 

4.1359e-4 

-0.18089 

1.0000 

0 

4 

4.6274e-4 

-0.18101 

0.3882 

1 

5 

2.4101e-4 

-0.18107 

0.2402 

1 

6 

2.6292e-4 

-0.18108 

0.0750 

2 

7 

8.9793e-5 

-0.18109 

0.0333 

3 

8 

2.3358e-6 

-0.18109 

0.0375 

2 


Table 4: Convergence behavior for Algorithm 3.11 Parameters: p = 4, g = 0.2 and 7 = 10 


the viscosity of shear-thickening materials increases with the rate of shear strain. In this case, since 
the shear stress transmitted by a fluid layer decreases toward the center of the pipe, the velocity 
takes a conical form with a flat part in the exact center of the geometry. 

In Table [4] we show the number of iterations that Algorithm |3. 11 needs to achieve convergence. 
We also show the evolution of \ J' fe (ufc)|/|J^, fe (lto)|, otk and the number of inner it¬ 

erations needed by Algorithm |4.3| to achieve convergence. The Algorithm performs as expected, 
i.e., the value of the functional is monotonically reduced at every iteration. Further, the residual 
behaves typically as in a deepest descent algorithm, but it shows fast local convergence. This be¬ 
haviour can be appreciated in Figure [7j This fact can be explained due to the stronger regularity 
of the differential operator when p > 2. As soon as g increases, this effect will be lost. This will be 
shown in the next experiment. 

Next, note that the step a *. does not have a monotone evolution during all the iterations. Also, 
the line search Algorithm |4.3j for some iterations, needs no inner iterations to achieve convergence. 
These facts can be explained due to the stronger convexity that the functional exhibits when p > 2, 
which implies that <Pk(ot) = J(uk + awk ) is better approximated by the quadratic model m*,. 


4.4.2 Experiment 2 

In this experiment, we set II to be the unit ball, and we compute the flow of a Herschel-Bulkley 
material with p = 10. We analyze the behavior of the algorithm with / = 1 and compare the 
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Figure 9: Velocity profiles for p = 10. 


performance of the Algorithm for g = 0.1 and g = 0.4. 

In Figure [8] the calculated velocities for <7 = 0.1 and g = 0.4 are depicted. As stated in the 
previous experiment, small values of g make the problem be close to the classical p -Laplacian 
problem. In this case, the Algorithm exhibits good performance. On the other hand, bigger values 
of g make the problem less regular. Also, from the mechanical point of view if the values of g 
increase, the size of the inactive zones increases as well. Therefore, the problem is more difficult to 
be approximated (see mm)- 

Regarding the performance of the Algorithm, in Figure [TO] the evolution of the error for g = 0.1 
and g = 0.4 are depicted. For g = 0.1, the error evolves in a typical way and the Algorithm achieves 
convergence in 14 iterations. On the other hand, for g = 0.4 the error is very oscillating and the 
Algorithm needs 27 iterations to achieve convergence. As expected, the fact that p and g increase 
provokes instabilities in the algorithm. 

Finally, in Table [5] we show the behaviour of the algorithm 3.11 in different meshes, considering 
a fixed value of g = 0.1. Here it can be appreciated that the Algorithm requires more iterations 
to achieve convergence as the mesh gets finer. This fact suggests that the Algorithm is not mesh 
independent. This is not shocking news, since the convergence result for Algorithm 3.11 


was 


obtained in a finite dimensional space. However, the algorithm still requires relatively few iterations 
to produce reliable solutions with low computational cost. 
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Figure 10: Calculated residual \ J' h (uk)\/\Ji y h(uo)\ for p = 10 and g = 0.1 (left) and g = 0.4 (right). 
Parameters: 7 = 10 3 . 


5 = 0.1 

hi 

h 2 

h 3 

/14 

Iter. num. 

l7,ftA)l/l7,»(^o)i 

J'y 

13 

5.242e-7 

-0.550987 

14 

3.868e-4 

-0.582629 

22 

3.931e-6 

-0.590403 

37 

4.518e-6 

-0.592738 


Table 5: Convergence behavior for Algorithm 3.5 Parameters: p = 10, g = 0.1 and 7 = 10 3 


4.4.3 Experiment 3 

One key issue in our approach is the size of the regularization parameter. In fact, theoretically, we 
obtain a better approximation for the problem when 7 is big. However, it is not a good strategy 
to directly run the Algorithms with high values for the parameter, since instabilities can arise in 
the process. In order to help the regularization parameter reach high values, we perform a simple 
but effective continuation technique: given 7 ^, we run the algorithm and obtain the corresponding 
solution • Next, we set 7^+1 = 107 &, initialize the algorithm with and run it to obtain 
it ^ . We stop this process when 7 equals 10 6 . 

As stated before, a challenging problem when using the p-Laplacian operator arises when p is 



Figure 11: Calculated velocity u for p = 100 and g = 0.3 at 7 = 10 5 . 
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Ik 

10 1 

10 2 

10 3 

10 4 

10 b 

10 B 

Iter. num. 

35 

11 

5 

9 

2 

1 


-0.2452 

-0.2354 

-0.2344 

-0.2343 

-0.2343 

-0.2343 

inbju 

9.0090 

9.0066 

9.0063 

9.0063 

9.0063 

9.0063 


Table 6 : Behavior of the continuation technique. Parameters: p = 100 and g = 0.4 


big. Therefore, we are interested in the computation of the flow of a Herschel-Bulkley material with 
p = 100 and g = 0.3. If we set 7 = 10 6 and run the Algorithm 3.11 convergence is not achieved. 
However, by using the continuation strategy, we obtain the solution for this problem. 

The convergence history is shown in Table [ 6 j We show the number of iterations that the 


Algorithm 3.11 needs to achieve convergence for each 7 ^, the value of the functional ~ity k ) and 


the norm of the calculated velocity IT 7 . It is possible to observe that, although the continuation 
technique helps the algorithm achieve convergence for high values of 7 the value of the functional 
and the norm of the velocity stabilize as soon as 7 *, equals 10 3 . This fact can be explained since 
the regularization procedure is sharp. Thus, for values around 7 ^ = 10 3 provides reliable results 
for the problem. However, we think that a further research in path following methods can clarify 
these aspects (see 


5 Conclusions 

In this paper, we focused on the numerical resolution of a class of variational inequalities of the 
second kind involving the p-Laplacian operator and the L 1 -norm of the gradient. The non differ¬ 
entiability of the associated functional was overcame with a Huber regularization procedure. This 
kind of local regularization has proved to be efficient in the context of this kind of problems. Based 
on optimization and variational techniques, we proposed preconditioned descent algorithms for 
coping the two cases 1 < p < 2 and p > 2. For the first case, we proposed an infinite dimensional 
descent algorithm and proved a global convergence result for it. The second case posed a difficult 
analytical issue, due to the lack of regularity of the candidates for descent directions. Thus, we 
proposed an algorithm in a finite dimensional setting and proved a global convergence result for 
this algorithm as well. Several numerical experiments were carried out to show the main features of 
the numerical approach. These numerical examples were constructed focusing on the applications 
to the flow of Herschel-Bulkley materials. Due to the structure of all the algorithms proposed, it 
was only necessary to solve one linear system at each iteration of the algorithms. This fact implied 
a low computational cost for all our numerical realisation. 

In order to continue this research, we consider that a deeper analysis of the case p > 2 is an 
interesting perspective. Here, the use of H s spaces provides a promising way to follow. Also, 
the combination of this approach with multigrid algorithms will be useful in order to cope more 
challenging problems, such as thep-Stokes problem (2D and 3D flows of Herschel-Bulkley materials). 
Finally, the analysis and simulation of blood flow models involving the Herschel-Bulkley structure 
looks like a very promising field of research. 
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